program harmonic
!--harmonic oscillator with Runge-Kutta method
use cdata
implicit  none
real :: y,t,f
end program

module cdata
implicit none
integer :: ndim
end module
!--------------------------------------------
subroutine eval_rhs(x,time,dx)
use cdata
implicit none
real,dimension(ndim) :: x,dx
real :: time
real,parameter :: omega=1.
intent,(in) :: x
intent,(out) :: dx
dx(1)=x(2)
dx(2)=-(omega**2)*x(1)
end subroutine
